Effect of the boundary condition on the vortex patterns in mesoscopic 
three-dimensional superconductors - disk and sphere 
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The vortex state of mesoscopic three-dimensional superconductors is determined using a minimiza- 
tion procedure of the Ginzburg- Landau free energy. We obtain the vortex pattern for a mesoscopic 
superconducting sphere and find that vortex hues are naturally bent and are closest to each other at 
the equatorial plane. For a superconducting disk with finite height, and under an applied magnetic 
field perpendicular to its major surface, we find that our method gives results consistent with pre- 
vious calculations. The matching fields, the magnetization and i/cs, are obtained for models that 
differ according to their boundary properties. A change of the Ginzburg-Landau parameters near 
the surface can substantially enhance Hc3 as shown here. 
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I. INTRODUCTION 

In the last decade the response of a mesoscopic su- 
perconducting disk to a perpendicular magnetic field has 
been theoreticall}*^!^ and experimentally^ studied. The 
small volume to surface ratio of mesoscopic supercon- 
ductors brings new and interesting physical properties 
such as giant vortices, recently detected thanks to new 
advances in small-tunnel-junction technology^. Previous 
studies of mesoscopic systems were based on two dimen- 
sional (2D) theory, where the superconducting conden- 
sate was assumed not to vary along the direction of the 
magnetic field. This assumption is not taken here and 
we study finite size extreme type-II mesoscopic supercon- 
ductors using a truly three-dimensional (3D) theoretical 
approach previously applied to a bulk superconductor-^'^ 



Only a few vortices fit inside a mesoscopic supercon- 
ductor whereas for a bulk superconductor, with non- 
superconducting inclusions inside, the number of vortices 
is uncountable. By inclusion we refer to a pinning center 
with the size equal to a few multiples of the coherence 
length, ^. These two systems have a similar properties 
because of their mesoscopic scale structure. For instance, 
giant vortices are naturally found in mesoscopic super- 
conductors but not in bulk superconductors, where the 
nucleation of a vortex line with multiple magnetic flux 
N^o is energetically forbidden and only the nucleation of 
N individual vortex lines with $o is possible (^o is the 
quantum of flux). However this picture does not hold in 
the presence of inclusions. For instance for a bulk super- 
conductor, Mkrtchyan and Shmidt^, Buzdin^, and some 
of us^'^^ have shown that a columnar defect can hold a 
multiple magnetic flux N^o- 

Metastability, matching fields, occupation numbers 
and giant vortices have been experimentally studied in 
2D bulk superconductors with inclusions, namely, su- 
perconducting films with an array of two-dimensional 
mesoscopic pinning centers consisting of not fully perfo- 



rated holes (blind holes) fully perforated holes (open 
holes) ^^^^^ and micro holes^^^^^. A similar 3D bulk super- 
conductor with a truly three-dimensional arrangement of 
internal inclusions has yet to be experimentally realized 
though it has been theoretically studied^' Such inter- 



nal inclusions are present as a random array 
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LREBaCuO superconductors (where LRE is a light rare 
earth element such as Nd, Sm, Eu and Gd) and some of 
the present ideas may be useful to explain their unusual 
properties. Such internal inclusions bring new features 
to vortex physics as just a single one can trap many vor- 
tex lines in its neighborhood. The regular array theo- 
retical study of the 3D bulk superconductor with inclu- 
sions was done in the context of a modified version of the 
Ginzburg-Landau (GL) theory. Other studies based on 
the Ginzburg-Landau theory for 3D systems have been 
done, including shells^^ and constricted superconducting 
wires^^ both in the extreme type-II limit. 

In this paper we apply the same theoretical approach 
to study a mesoscopic superconductor. A 3D lattice of 
inclusions turns into a 3D lattice of mesoscopic super- 
conductors with the same geometry when the insulat- 
ing regions are replaced by superconducting ones and 
vice-versa. The 3D lattice of mesoscopic superconduc- 
tors becomes a set of individually equivalent mesoscopic 
superconductors for a London penetration length much 
larger than their size. In this case the local field is con- 
stant and equal to the applied field everywhere. In this 
way we obtain the vortex patterns for a single meso- 
scopic superconductor, here of a sphere and of a disk. In 
the past the vortex patterns of mesoscopic superconduc- 
tors were obtainedi^^i^UiSi in the limit of an extremely 
thin film. This condition renders the variation of the 
Cooper pair density along the magnetic field negligible. 
Baelus and Peeters^l studied several different flat ge- 
ometries typically with thickness 0.1^ and obtained the 
vortex patterns from the two-dimensional GL equations 
supplemented by the Saint-James-de Gennes^^ bound- 
ary conditions at the edge. They considered a Ginzburg- 
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Landau parameter, the ratio of the London penetration 
to the coherence length, = 0.28, and solved the two 
GL equations. Here we study a thick disk and compare 
our results to theirs^^ taking into account that in their 
case the magnetic response to an external applied field 
is much stronger than here. We only report results for 
Hi ^ oo although our method is not restricted to this 
limit. Notice that for our case the boundary conditions 
are truly three-dimensional, and so, imposed in all direc- 
tions including perpendicularly to the fiat geometry. 

The major new results of this paper can be summarized 
as follows, (i) We find the vortex pattern for a meso- 
scopic sphere, with radius Rg = 4.0^, a problem whose 
solution is beyond the scope of previous 2D techniques, 
(ii) We show that a slight change of the Ginzburg-Landau 
parameters near the edge can substantially increase the 
Hc3 field. A thin layer covers the superconductor and 
separates it from the outside insulating world. This layer 
is also superconducting but with effective GL parameter 
slightly different from those inside. For a bulk system 
the phenomenological GL parameters are known to be 
related to the microscopic parameters in the following 
way: ao ~ {kTcY /ep^ /3 ~ (/cTc)^/(ei?^), where Tc is the 
critical temperature, ei? is the Fermi energy, and n is the 
electronic density. Similar relations should exist in case 
of a mesoscopic superconductor although we don't obtain 
them here. We just show that a decrease near the edges 
of the effective Cooper pair mass, m, and of ao, lead to 
an enhancement of Hcs. Therefore the present approach 
is interesting for a system with a small volume to sur- 
face ratio because there a slight change at the boundary 
over a distance less or equal to ^ is found here to make 
a significant difference. The present approach relies on a 
free energy minimization procedure carried in the whole 
space, including the world outside the superconductor, 
where the order parameter is found to vanish. The de- 
cay of the Cooper pair density at the boundary, from a 
finite value inside the mesoscopic superconductor to zero 
outside, is treated here. Notice that standard differen- 
tial equation approaches, such as that of Ref. [2l|, only 
treat the volume internal to the superconductor, and do 
not treat the order parameter discontinuity at the edge, 
from a finite value to zero at the outside world. In this 
paper we study three kinds of boundary conditions and 
discuss them in the context of a disk of radius R = 4.0^. 

Below we provide a short description of the disk and 
sphere boundary problems treated here. We chose to give 
them names that recall their major features: (i) sharp: 
a disk is considered and its boundary treatment is the 
standard one used for comparison with all other models. 
A coarse grained grid is used and gives a fast and effi- 
cient convergence to the final configuration. The vortex 
states are satisfactorily described here. Its name stems 
from the sharp definition of the edge, (ii) mesh: this 
model is the same as sharp except with a refined grid, 
which contains 8.2 times more grid points, (iii) sphere: 
a sphere is treated here with the same grid coarseness 
along the disk radial direction as in the sharp model, (iv) 



BP2D: this is the disk reported by Baelus and Peeters21 
using their two-dimensional approach, (v) smooth: this 
type of boundary was previously used in Refs. [23ll23 for 
insulating pinning spheres inside a bulk superconductor 
and is used here for the disk. Its major property is that 
the supercurrent normal to the surface does not disap- 
pear abruptly but over some small region (fraction of ^) . 
(vi) step: This model contains a superconducting layer 
that sets the disappearance of superconductivity. Thus 
there are two concentric disks and we find that this in- 
termediate layer stabilizes the superconducting state in 
the inner disk. This model features a very high i^cs as 
compared to the other models. 

The paper is organized as follows. In Section [III we de- 
scribe our theoretical approach valid for the following two 
complementary situations: (i)superconductor with non- 
superconducting inclusions and (ii) mesoscopic supercon- 
ductor. In Section IIIII we describe the disk and sphere 
boundary models and discuss their properties obtained 
through our numerical simulations. In Section [iVl we 
compare the models and discuss many of their common 
features. Finally in Section |Vl we summarize the main 
achievements of this work. 



II. THEORETICAL ASPECTS 

One of the advantages of the present method is that it 
can easily incorporate the shape of the mesoscopic super- 
conductor, which is done at the free energy level, given 
by the density expansion below. 



F 
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where D = {h/i)V — qA/c^ q is the Cooper pair charge, 
and r(r) is a step-like function, equal to one inside the 
mesoscopic superconductor and zero outside. The r(r) 
contains the geometry of the mesoscopic superconductor. 
The corresponding Euler-Lagrange equations are given 
by, 



f-Vr(r) • +r(r) 



2m 



^T{r)ao(T - T,)^ + = 0, (2) 



\/xh= ,J=J^r r 

c 2m 



i/j'^DiP^iDipyilj .(3) 



These modified GL equations automatically incorporate 
the appropriate boundary conditions through the step- 
like function T(r), discontinuous at the edge, equal to one 
inside, and zero outside, for the mesoscopic superconduc- 
tor. The gradient of r(r) is zero everywhere except at the 
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interface, where it diverges. Any finite and physical so- 
lution must obey Vr • D^Ij = because this divergence 
is proportional to the Dirac delta function. For instance, 
along the radial direction of the disk: r(r) = 1 for r < 
and r(r) = for x > thus the derivative becomes 
Vr = rdr{r)/dr = —fS{r — R). Let /(r) be any function 
describing products of the order parameter and its deriva- 
tives. The product f{r)dr{r)/dr diverges at the border 
and the only way to make it vanish there is through the 
condition f{R) = 0. Thus the well-known Saint- James- 
de Gennes^^ boundary condition, n ■ DiIj\^ = 0, is re- 
covered here. As the superconducting order parameter 
is defined everywhere in the unit cell, including outside 
the mesoscopic superconductor, where Eq. ([1]) becomes 
F = (1/^) / M (3\iIj\^/2 + P/Stt ] outside the super- 
conducting volume. It must vanish as a result of the 
free energy minimization, and variation with respect to 
and A shows that the minimum is reached for = 
and V X h = according to Eqs. ([2])-(|3]). It is possi- 
ble to obtain more elaborate versions of the GL theory, 
such as the one containing a local depression of the criti- 
cal temperature through a function Tc(r)^ in Eq. ([2|). In 
this work the free energy is normalized by the constant 
Fq = H^/Sir and all fields are normalized in terms of 
the upper critical field, i^c2- Lengths are in units of the 

coherence length, ^(T) = ^J /2mao{Tc — T), and the 

density is normalized by {ao{Tc — T)//3)^ such that 
its maximum value of 1 is reached, for instance, inside 
a bulk superconductor (no boundaries) for zero applied 
field. 

We stress some differences in the application of the 
present method to the two complementary problems. For 
the former case the magnetic induction, B = J dv h/V, 
is constant whereas for the latter the applied field H is 
constant. In the former case the unit cell edges are fully 
inside the superconductor and this introduces into the 
theory integers associated to the periodic boundary con- 
ditions imposed by the unit cell. These integers follow 
the condition that the order parameter be single- valued. 
Though IV^p and h are periodic, ip and A only need to 
coincide at the unit cell surfaces up to a gauge transfor- 
mation, whose expression gives room to introduce these 
integers. 

V^(r + L^) = e^^^^^V(0 (4) 
A(r + 4) = A(r)+VA^(f) (5) 

where is the unit cell length, A^(r) is a scalar gauge 
function and ji = y oi z. The minimization procedure 



shows that such integers are nothing but the number of 
vortices in the unit cell, and the magnetic induction is 
fully determined by them. However for the complemen- 
tary problem, the mesoscopic superconductor is fully in- 
side the unit cell and its boundaries are away from the 
unit cell edges. Consequently there is no single-valued 
condition on the order parameter and so, these integers 
do not exist at all. Consequently, the independent ther- 
modynamic field in this case, which is the applied field 
varies continuously. 

Since there are no screening currents the local field, de- 
fined as = V X A, is the external applied field h = H. In 
this large limit the magnetization is directly determined 
from M = const J dv f x where J is the supercur- 
rent. An extra condition determines the remaining free 
parameter const ^ and consequently the demagnetization 
constant D of the mesoscopic superconductor: for small 
J?, that is, in the Meissner phase, we impose the condi- 
tion that H ^AitDM = 0. In contrast, in the approach of 
Baelus and Peeters^^ for finite the magnetization is di- 
rectly obtained from the difference between the magnetic 
induction and the applied field. 

The minimization of the GL free energy, done nu- 
merically through the so-called Simulated Annealing 
method^i^^, is carried in a discrete three-dimensional 
space. The free energy is adapted to keep its gauge 
invariance in this discrete space. A cell, that consists 
of an orthorhombic box containing N^.Ny.Nz points for 
this purpose. Every point in this cell, belonging to 
the mesoscopic superconductor or not, has associated 
to it the fields V^(na;, n^y, n^), and A(naj, n^y, n^), where 
n^ l,...,A^a,, ny l,...,Ny, and n^ 1,...,7V2. 
The physical volume of the box is ■ Ly ■ Lz^ where 
Lx Gix{Nx - 1), Ly ay{Ny - 1) and a^^N^ - 1). 
The distance between two consecutive points along the 
axes of the box is a^^, a^, and a^. The discrete theory, 
given by Eq. (|6]), properly describes the properties of the 
continuous theory under the condition that ^ be much 
larger than a^^, a^y, and a^, the grid resolution. 

In the discrete free energy, given by Eq. ([6]), grid 
points outside and inside the mesoscopic superconduc- 
tor are coupled through gradient terms. For instance in 
case of no applied field, this coupling is proportional to 
Wout ~ i^in\^' The fact that ip'^^^ causes ipin to get a 
lower value than deep inside the sample, where the kinetic 
energy vanishes as the order parameter is the same in all 
points. For this reason the density never reaches its 
maximum bulk value due to the small volume to surface 
ratio. 
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III. COMPARISON OF THE DIFFERENT 
MODELS 

The models introduced in Section [H have their prop- 
erties summarized in Tables |T] and |TT1 and some of their 
free energy and magnetization properties are described 
in Tables [III] and HVl respectively. The function r(r) 
is taken as the product of independent orthogonal direc- 
tion functions for all models studied here. Below is a 
summary of their features, (i) sharp: This model treats 
the boundaries of a disk of radius R = 4.0^ and height 
d = 1.0^ through a r(r) function, 1 inside and out- 
side, both along the radial and the major axis direction: 
r(r) = Tp{p)-Td{z). (ii) mesh: The same disk and bound- 
aries of the sharp model is treated here but with a denser 
grid, 61 X 61 X 26 instead of 41 x 41 x 7. (iii) sphere - 
This model treats a sphere with stepwise r(r) function 
such as in the sharp model, (iv) BP2D: This is the disk 
of Ref. [m. It has a 0.1^ height and a 128 x 128 x 1 
grid is used. Although of its two-dimensional treatment 
it contains 1.4 times more grid points than the three- 
dimensional sharp disk model, (v) smooth: For this 
model a larger height is taken, 2.0^, to help stabilize 
the order parameter inside the disk. The smoothness 
of Td{z)^ which reaches values below 1 inside the region 
\z\ < (i/2), tends to downgrade the order parameter in- 
side the disk. For our numerical study we selected, for the 
exponential parameter of Table [U N=8. (vi) step: The 
height is d = 1.5^ and the r(r) function varies stepwise, 
taking values 0, 0.8 and 1. The choice of intermediate 
value 0.8 is rather arbitrary and we have found that low- 
ering this intermediate value to 0.5, for instance, causes a 
substantial increase of i^c3, as compared to here. Thus a 
drop of r near the border, and so of the corresponding GL 
parameters, can severely affect i^cs- Fig. [1] shows the 
normalized density for zero applied field (Meissner phase) 
versus the distance from the geometric center of the disk 
along the radial direction, and in case of the sphere, this 
distance is along the radial direction inside the equato- 
rial plane. For clarity the six models considered here were 
split into two subsets shown in different plots. For com- 
parison purposes the BP2D model is shown in both plots 
(red). Notice that for the BP2D model, as well as for the 
sphere, the maximum density is 1.0, but not for the other 
disk models whose maximum density is about 0.8. The 



sphere has a larger volume than the disk, and so, the 
surface is not so effective to alter the order parameter in 
its center. 

In presence of an applied field the numerical simulation 
is carried in the following way. For zero applied field a 
random configuration of the order parameter is assumed 
inside the cell and a search for the minimum of the free 
energy is carried out. The applied field is increased at 
constant steps and for a given field one assumes as the 
starting order parameter configuration the one found for 
the previous field. This procedure is carried sequentially 
until the last critical field, i^c3, is reached and the order 
parameter vanishes everywhere. Next the applied field is 
lowered backwards to zero applied field. A typical feature 
of mesoscopic superconductors'^^ is a saw-tooth structure 
for the descending field of magnetization curve. The two 
curves do not coincide, the ascending one has a stronger 
diamagnetic signal than the descending curve. The meso- 
scopic superconductor exhibits hysteresis as observed, e. 
g., in Al disks^, and theoretically obtained in previous 
studies''^'. Metastability is also found in the complemen- 
tary problem of bulk superconductors with inclusions^i^. 

Notice that all the magnetization curves shown 
here (Figs. [21 O [4j) decompose into independent non- 
intersecting lines and result from ancestor curves that 
contain their ascending and descending branches. The 



TABLE I: The different models considered in this paper where 
r(r) = Tp{p) • Td{z). Rs is the sphere radius. Ri and di are 
the internal disk radius and height, respectively. 



Model 






sharp 




1 P<R 
p> R 


-\ 


1 2\z\<d 
2\z\ > d 


mesh 


idem 


idem 


sphere 


r(f) = 


[ 1 r<Rs 
[0 r>Rs 




BP2D 


dif. eq. 




smooth 




Td = 2 




step 


...| 


p<R^ 
0.8 Ri<p<R 
^0 p>R 


■H 


^1 2\z\<d^ 

0.8 di<2\z\<d 
2\z\>d 



5 



TABLE II: The parameters of the different models used in 
our numerical calculation. 



model 


grid"" 


cell^ 


parameters^ 


sharp 


(41,41,7) 


(0.3,0.3,0.5) 


d=1.0 


mesh 


(61,61,26) 


(0.2,0.2,0.2) 


d=1.0 


sphere 


(41,41,41) 


(0.3,0.3,0.3) 


Rs=4.0 


BP2D 


(128,128,1) 


(0,0,-) 




smooth 


(41,41,13) 


(0.3,0.3,0.5) 


d=2.0,N=30 


step 


(41,41,7) 


(0.3,0.3,0.5) 


d=1.5,R^=3.5, di=0.5 



^The number of grid points for the three cell directions: 

^The lattice spacing for the three cell directions: (a^, ay,az). 
^All the lengths are in units of ^. The disk radius is R = 4.0 for 
all cases. 



ascending field the vortex pattern moves from L inde- 
pendent vortices at low field to giant vortex states at 
high field, whose total angular momentum must add to 
L. Therefore the present method is able to reproduce 
the well-known features of mesoscopic superconductors 
found by other authors using different appro achesi*^. 

Table [TTTl shows the matching fields Hll+i between 
two nearest angular momentum states for the six models 
considered here. Table HVl is useful for model comparison, 
as it shows the maximum value of —A-kDMl/ Hc2 for each 
L lines and its corresponding applied field Hl/ Hc2 • 



IV. DISCUSSION 



saw-tooth structure is a sum of segments, which are parts 
of the independent non-intersecting lines. Two consecu- 
tive segments are connected by a vertical jump. The rea- 
son for such vertical jumps resides in the free energy curve 
which also consists of independent but intersecting lines. 
In fact these intersections define the so-called matching 
fields which correspond to cross-sections of free energy 
lines of neighboring states. Above the matching field 
the free energy of the higher state is lower than that of 
the preceding one and thus the higher state is preferred. 
This is also the reason for the saw-tooth character of the 
hysteresis curve. Depending on how the numerical pro- 
cedure is carried (the magnetic field step, the simulated 
annealing temperature, etc) one obtains a different saw- 
tooth structure that always falls over the same set of inde- 
pendent non-intersecting lines. The presence of distinct 
lines in the magnetization and free energy curves reveals 
a parameter that remains constant upon field sweep. A 
look at the order parameter phase reveals that it takes 
variations from to 2tt and the number of such varia- 
tions remains constant throughout a magnetization line. 
Thus this parameter is the total angular momentum L^'^ . 
Along any of these lines the angular momentum remains 
constant, such that each line can be labeled by L. For 



TABLE III: The matching fields h^L+i between the angular 
momentum states L and L+1 for the different models consid- 
ered here. 



hii+i 


sharp 


mesh 


sphere 


BP2D 


smooth 


step 


hoi 


0.31 


0.30 


0.41 


0.39 


0.30 


0.31 


hi 2 


0.51 


0.51 


0.65 


0.59 


0.52 


0.54 


h2 3 


0.69 


0.68 


0.84 


0.74 


0.69 


0.72 


h3 4 


0.85 


0.83 


1.00 


0.89 


0.84 


0.88 


h4 5 


0.99 


0.98 


1.15 


1.02 


0.99 


1.03 


hse 


1.14 


1.12 


1.28 


1.16 


1.13 


1.20 


he? 


1.28 


1.26 


1.42 


1.30 


1.27 


1.34 


h7 8 


1.43 


1.41 


1.56 


1.43 


1.41 


1.48 


hsQ 


1.57 


1.54 




1.57 


1.54 


1.62 


hg 10 


1.70 


1.69 




1.71 


1.68 


1.84 


hio 11 


1.87 


1.82 




1.84 




2.05 


hii 12 












2.19 



In this section we compare all models to the sharp 
boundary model for the mesoscopic disk whose proper- 
ties are shown in Tables IIIII and llVi The effect of the 
number of grid points in our calculations can be checked 
in Fig. [2] as the mesh model has 8.2 times more grid 
points than the sharp model. The mesh model has a 
lower free energy and a higher magnetization than the 
sharp model, but their values differ by less than one per 
cent. Effects due to the grid become only noticeable for 
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FIG. 1: (Color online) Cooper pair density |^|^ vs. the dis- 
tance from the center of the disk for the case of the Meissner 
state in the absence of an applied field. The symbols corre- 
spond to 1^1^ values at the mesh grid points. 
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TABLE IV: Points {Hl/ Hc2^ —AtiDMl / Hc2) are the maximum of the magnetization hnes associated to the angular momentum 
L curves for the different models. 



L 


sharp 


mesh 


sphere 


BP2D 


smooth 


step 





(0.38,0.25) 


(0.39,0.25) 


(0.43,0.28) 


(0.58,0.44) 


(0.39,0.25) 


(0.44,0.28) 


1 


(0.58,0.19) 


(0.59,0.20) 


(0.66,0.19) 


(0.73,0.39) 


(0.59,0.19) 


(0.64,0.22) 


2 


(0.75,0.14) 


(0.75,0.15) 


(0.84,0.13) 


(0.86,0.34) 


(0.75,0.14) 


(0.81,0.17) 


3 


(0.91,0.11) 


(0.90,0.12) 


(0.98,0.090) 


(0.98,0.29) 


(0.90,0.10) 


(0.98,0.13) 


4 


(1.05,0.082) 


(1.04,0.089) 


(1.14,0.060) 


(1.11,0.24) 


(1.04,0.077) 


(1.13,0.11) 


5 


(1.19,0.060) 


(1.18,0.067) 


(1.27,0.038) 


(1.23,0.20) 


(1.17,0.055) 


(1.27,0.083) 


6 


(1.32,0.044) 


(1.30,0.049) 


(1.39,0.022) 


(1.34,0.15) 


(1.30,0.038) 


(1.41,0.064) 


7 


(1.44,0.030) 


(1.43,0.034) 


(1.50,0.011) 


(1.45,0.12) 


(1.41,0.025) 


(1.55,0.048) 


8 


(1.56,0.019) 


(1.54,0.022) 


(1.60,0.0034) 


(1.56,0.083) 


(1.53,0.015) 


(1.67,0.036) 


9 


(1.67,0.011) 


(1.67,0.013) 




(1.66,0.054) 


(1.64,0.0072) 


(1.80,0.026) 


10 


(1.78,0.0046) 


(1.77,0.0057) 




(1.78,0.029) 


(1.74,0.0023) 


(1.92,0.016) 


11 


(1.88,8.8x10"^) 


(1.86,9.7x10"^) 




(1.87,0.0096) 
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intermediate fields, but not in the single vortex region. 
The comparison between sharp and mesh shows that the 
present numerical approach is robust and displays very 
little quantitative dependence on the grid. Fig. [2] also 
shows the iso-density three-dimensional plots of four typ- 
ical vortex configurations selected to display 1, 2, 3 and 4 
vortex states, respectively. Their corresponding applied 
field, magnetization and free energy values can be read 




FIG. 2: (Color online) The sharp (red) and the mesh disks 
free-energy and magnetization curves are shown here. Iso- 
density plots for selected applied fields are shown here to illus- 
trate the first four cases of vortex patterns found for the thick 
disk. The three-dimensional figures are iso-contours taken at 
20% of the maximum density IV^I^. Each iso-contour is a single 
surface, the sum of the vortices and the external surface. 



from Fig. O 

Fig. [3] shows results of the magnetization and free 
energy for the sphere model, and also for the sharp 
disk model. Three-dimensional iso-density and two- 
dimensional contour plots are shown for four selected 
vortex configurations, whose location is indicated in the 
free-energy curve. These four vortex configurations illus- 
trate general features of the vortex lines inside the sphere. 
A vortex line must reach the surface perpendicularly in 
order to avoid a supercurrent component pointing out- 
wards the surface^^^^^. Because of the small volume to 
surface ratio of this Rs = 4.0^ sphere, the vortex lines 
are strongly affected by this surface effect and as a result 
they are curved everywhere inside the sphere. The lines 
are closely packed in the equatorial plane as shown by 
the three-dimensional iso-density plots. These plots are 
drawn at 20% of the maximum order density and each 
plot consists of a single iso-density surface. The north 
pole part of these plots provide a view of the vortex be- 
havior at the surface, but the translucent properties of 
these three-dimensional plots make it difficult to have 
the same view at the south-pole. Fig. [3] also shows two- 
dimensional contour plots associated to two selected cuts 
of the sphere, taken at the equatorial (half-plane) and 
half-way between the north (or south) pole and the equa- 
tor planes (quarter-plane). These contour plots contain 
ten contour regions, shown in different colors, ranging 
from maximum density (red) to minimum density (blue) . 
They also show that the vortex lines are closely packed 
at the equatorial plane and also that the vortex core is 
larger near the surface than inside the sphere. The sphere 
has a stronger magnetic signal as compared to the disk 
for low fields, but for high fields up to i^cs the situation 
turns and the disk acquires a stronger signal. In fact the 
sphere only supports 9 vortex states whereas the disk 12 
states. As the field increases the vortex configuration 
in case of the sphere disappears faster than in the disk, 
probably due to the existence of vortex lines of different 
lengths. 
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FIG. 3: (Color online) The sphere and the sharp (red) disk free-energy and magnetization curves are shown in the right panels. 
Three-dimensional iso-density plots and two dimensional density contour plots of IV^I^ are also shown in the right panels. Each 
three-dimensional iso-contour is a single surface, made of the sum of the vortices and the external surface. Two-dimensional 
contour plots are taken at the half (equatorial) and at the quarter plane that cuts the sphere perpendicularly to the applied 
field direction. 



Fig. m shows comparative analysis of the free energy 
and magnetization curves to the sharp model for the 
BP2D^ step^ and smooth models, {i) BP 2D -sharp: Their 
different K: values yield significantly different magnetiza- 
tion curves. The BP 2D model has a very strong dia- 
magnetic signal lower free energy states, because it has 
a more effective shielding to the applied field. Neverthe- 
less the models show qualitative similarities. They both 
have the same number of 12 angular momentum states, 
as shown in Table IIIH and present a fair agreement be- 
tween matching fields in the high field region. This is 
explained by the weakening of the diamagnetic currents 
for high field that turns the {BP 2D) similarly to the sharp 
model, {ii) step- sharp: The presence of an intermediate re- 
gion at the boundary enlarges surface effects as compared 
to the sharp case. The diamagnetic response is stronger, 
and the free energy is lower, in all L lines and in fact it 
allows for two extra vortex states, according to Table Hill 



These features are not a consequence of a slight differ- 
ence in height between the two models. {iii) smooth- sharp: 
The smooth model treats the boundary differently from 
sharp in case the smooth r function decay takes place 
over a distance larger than the mesh parameters a^^, a^, 
and a^. This is the case here but we find no substan- 
tial change in behavior by using the smooth model. This 
boundary was extensively used in previous problems of a 
superconductor with a periodic array of inclusions^^^. 

V. CONCLUSION 

The vortex patterns of truly three-dimensional meso- 
scopic superconductors, namely a disk and a sphere, were 
analyzed. They were obtained by numerical minimiza- 
tion (Simulated Annealing) of a modified GL free energy 
that already incorporates the boundary conditions. This 
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FIG. 4: a) and b): The two-dimensional BP2D disk and the three-dimensional s/iarp (red) disk magnetization and free-energy, 
c) and d): The step and the sharp (red) disk magnetization and free-energy, e) and f): The smooth and the sharp (red) disk 
magnetization and free-energy. 



procedure provides an efficient way to obtain vortex pat- 
terns in mesoscopic superconductors and needs relatively 
few grid points. The method is stable under changes 
of the grid size, and for a two-dimensional disk it repro- 
duces results of disk geometry previously studied by other 
methods^^). We find that slight changes of the boundary 
conditions, like the creation of a surface layer, increases 
the upper critical field and allows for an increase in the 
number of angular momentum states. In case of a meso- 
scopic sphere we find that the vortex lines are naturally 
curved due to strong surface effects as was recently also 
found in a wire with a constrictioni^i. 
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